On this page

Exercise 5-2: Oscillations in a Vertical Plane on Three PlatformsΒΆ

Here particles oscillate through the center within a vertical plane. The period for each platform is displayed once it finishes.

Index 0 is simple harmonic motion - Index 1 is a paraboloid - Index 2 is a concave-up hemisphere

Code translated from GW-BASIC provided in Exercise 5.2 of Stommel and Moore (1989)

"""
author: Victoria McDonald
email: vmcd@atmos.washington.edu
website: https://github.com/torimcd/coriolis-sm

"""
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.animation import FuncAnimation
import numpy as np
from IPython.display import HTML
%matplotlib inline

Here we set up the constants and initialize the global variables

# set up global variables
dt = 0.01   # time step
t = 0
g = 1       # gravity
d = 1
k = 1       # spring constant
k2 = k**2
c = k2/2
a = k2/g
mu = 0    # angular momentum
r = {}
dr = {}
ddr = {}
p = {}
dp = {}
x = {}
y = {}
flag = {0:0, 1:0, 2:0}
drl = {0:-1, 1:-1, 2:-1}
zzz = {}

x_plane = []
y_plane = []
x_para = []
y_para = []
x_hemi = []
y_hemi = []
l_plane = []
l_para = []
l_hemi = []
m_plane = []
m_para = []
m_hemi = []


# set up the surfaces
for i in range(3):
    r[i] = 0.5
    dr[i] = 0
    p[i] = 0

Here we set up the figure that the particles will travel on

%%capture

fig, axes = plt.subplots(2,3)

# turn the axes off
axes[0,0].axis('off')
axes[0,1].axis('off')
axes[0,2].axis('off')
axes[1,0].axis('off')
axes[1,1].axis('off')
axes[1,2].axis('off')

# set the aspect ratio
axes[0,0].set_aspect('equal')
axes[0,1].set_aspect('equal')
axes[0,2].set_aspect('equal')

# create outline
outline_1 = patches.Ellipse((0,0), 1, 1, 0, linewidth=2, fill=False)
outline_2 = patches.Ellipse((0,0), 1, 1, 0, linewidth=2, fill=False)
outline_3 = patches.Ellipse((0,0), 1, 1, 0, linewidth=2, fill=False)

# add outline to axes
axes[0,0].add_patch(outline_1)
axes[0,1].add_patch(outline_2)
axes[0,2].add_patch(outline_3)

axes[0,0].set_xlim(-0.6, 0.6)
axes[0,0].set_ylim(-0.6, 0.6)

axes[0,1].set_xlim(-0.6, 0.6)
axes[0,1].set_ylim(-0.6, 0.6)

axes[0,2].set_xlim(-0.6, 0.6)
axes[0,2].set_ylim(-0.6, 0.6)

axes[1,0].set_xlim(-1, 1)
axes[1,0].set_ylim(-1, 1)
axes[1,1].set_xlim(-1, 1)
axes[1,1].set_ylim(-1, 1)
axes[1,2].set_xlim(-1, 1)
axes[1,2].set_ylim(-1, 1)

plt.suptitle("Three Surfaces of Revolution")
axes[1,0].text(-0.5, -0.5, "plane")
axes[1,1].text(-0.5, -0.5, "paraboloid")
axes[1,2].text(-0.5, -0.5, "hemisphere")

Now we set up the particle objects to be updated in each step of the animation

# set up the particles
plane = axes[0,0].plot([], [], color='green')[0]
paraboloid = axes[0,1].plot([], [], color='red')[0]
hemisphere = axes[0,2].plot([], [],  color='gold')[0]

plane_section = axes[1,0].plot([], [], color='green')[0]
paraboloid_section = axes[1,1].plot([], [], color='red')[0]
hemisphere_section = axes[1,2].plot([], [], color='gold')[0]

This is the update() function that calculates the new position for each surface, and updates the particle objects at each step

n=500
def update(frame):
    global dt
    global g
    global d
    global k
    global k2
    global c
    global a
    global mu
    global r
    global dr
    global ddr
    global p
    global dp
    global x
    global y
    global t

    global x_plane
    global y_plane
    global x_para
    global y_para
    global x_hemi
    global y_hemi
    global l_plane
    global l_para
    global l_hemi
    global m_plane
    global m_para
    global m_hemi
    global flag
    global drl

    # draw cross sections first
    if frame < 70:
        i_map = np.arange(-0.7, 0.7, 0.02)
        i = i_map[frame]
        
        l_plane.append(i)
        m_plane.append(0)

        l_para.append(i)
        m_para.append(c*(i**2))

        if (a**2 - i**2) <= 0:
            plane_section.set_data(l_plane[:], m_plane[:])
            paraboloid_section.set_data(l_para[:], m_para[:])
            hemisphere_section.set_data(l_hemi[:], m_hemi[:])

        else:
            l_hemi.append(i)

            m_hemi.append(a-np.sqrt(a**2 - i**2))
                
            plane_section.set_data(l_plane[:], m_plane[:])
            paraboloid_section.set_data(l_para[:], m_para[:])
            hemisphere_section.set_data(l_hemi[:], m_hemi[:])
            

    else: 

        # integrations for plan view
        ddr[0] = mu**2/(r[0]**3)-k**2*r[0]
    
        aaa = (1+(2*c*r[1])**2)
        ccc = mu**2/(r[1]**3)
        bbb = -2*g*c*r[1] - r[1]*(2*c*dr[1])**2

        ddr[1] = (bbb+ccc)/aaa

        aa = -g*r[2]/(np.sqrt(a**2 - r[2]**2))
        bb = mu**2/(r[2]**3)
        cc = -dr[2]**2*(a**2*r[2]/((a**2-r[2]**2)**2))
        dd = 1+r[2]**2/(a**2*r[2]**2)

        ddr[2] = (aa+bb+cc)/dd

        for i in range(3):
            dr[i] = ddr[i]*dt + dr[i]
            r[i] = dr[i]*dt + r[i]

            #for j in range(3):
            dp[i] = mu/(r[i]**2)
            p[i] = dp[i]*dt + p[i]

            #for k in range(3):
            x[i] = r[i]*np.cos(p[i])
            y[i] = r[i]*np.sin(p[i])

            #for l in range(3):
            if flag[i] == 1:
                drl[i] = dr[i]
            else:
                zzz = dr[i]*drl[i]
                if zzz<0:
                    flag[i] = 1
                    axes[1, i].set_title("time("+str(i+1)+"): " + str(2*t)[0:5], loc='left')
                else:
                    drl[i] = dr[i]

        t = t+dt

        x_plane.append(r[0]*np.cos(p[0]))
        y_plane.append(r[0]*np.sin(p[0]))

        x_para.append(r[1]*np.cos(p[1]))
        y_para.append(r[1]*np.sin(p[1]))

        x_hemi.append(r[2]*np.cos(p[2]))
        y_hemi.append(r[2]*np.sin(p[2]))

        plane.set_data(x_plane[:], y_plane[:])
        paraboloid.set_data(x_para[:], y_para[:])
        hemisphere.set_data(x_hemi[:], y_hemi[:])

Finally, we construct the animation and embed it in this page

# Construct the animation, using the update function as the animation director.
animation = FuncAnimation(fig, update, frames=n, interval=10, repeat=False, blit=False)

# convert to a video to be embedded in web page
HTML(animation.to_jshtml())